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Abstract 

Normal surfaces are a key tool in computational knot theory and 3-manifold topol- 
ogy, and have featured in significant computational breakthroughs in recent years. Despite 
this, there has been little practical progress on algorithms that use fundamental normal 
surfaces, which are described in terms of a Hilbert basis for a pointed rational cone on 
a high-dimensional integer lattice. In this paper we develop and implement several algo- 
rithms to enumerate fundamental normal surfaces, by merging domain-specific techniques 
from normal surface theory with classical Hilbert basis algorithms. The most successful 
of these combines a maximal admissible face decomposition with the primal Hilbert basis 
algorithm of Bruns, Ichim and Koch, and in many cases can solve 168-dimensional enu- 
meration problems (based on 24-tetrahedron knot complements) in a matter of hours. As 
an application, we use this new algorithm to compute 164 previously unknown crosscap 
numbers in the Knotlnfo database of knot invariants. 

Keywords Normal surfaces, fundamental surfaces, algorithms, 3-manifolds, knots, cross- 
cap number 

1 Introduction 

Normal surfaces are a vital tool for algorithms and complexity in 3-manifold topology and 
knot theory. For algorithms, they allow us to reduce difficult topological searches to discrete 
problems on integer points in rational cones [531 13H] • For complexity, they allow us to produce 
polynomial-sized certificates that prove that unknot recognition and 3-sphere recognition are in 

NP mils]. 

Normal surfaces were introduced by Kneser |36| , and later developed for algorithmic use by 
Haken [52]. The key ideas are as follows. We present the input to our topological problem as 
a 3-manifold triangulation 7" (so, if the input is a knot, we triangulate the knot complement 
[23]), and we frame our problem as a search for some embedded surface in T with a particular 
property (for instance, for the unknot recognition problem we search for a spanning disc in the 
knot complement). We then show that we can restrict this search to normal surfaces in T, which 
are properly embedded surfaces that intersect the tetrahedra of T in a well-behaved fashion. 

Normal surfaces have the useful property that we can encode them arithmetically. Specifi- 
cally, if T contains n tetrahedra, we can encode normal surfaces as integer points in the normal 
surface solution cone which is a pointed rational cone in M7^. For many topological prob- 
lems, we can construct a finite list ^ of points in this cone so that, if there is some surface in 
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T with the desired property, then there is one described by a point in our Hst I£ . A typical 
algorithm would then build the cone construct the finite list I£ ^ rebuild the corresponding 
surfaces, and test each for the desired property. 

A key point of difference between topological algorithms is how we construct this finite list 

(i) For some algorithms, ^ contains the smallest integer point on each extremal ray of the 
cone ^ (sometimes their doubles are also required) . The corresponding surfaces are called 
vertex normal surfaces. 

Topological algorithms of this type include unknot recognition [22l[23], 3-sphere recognition 
[3T| |42] , computing knot genus in homology 3-spheres [47] , connected sum decomposition 
[31] , and Hakenness testing [SO] |34] . 

(ii) For other algorithms, ^ is the Hilbert basis for the cone ^; that is, all integer points 
in that cannot be expressed as a non-trivial sum of other integer points in The 
corresponding surfaces are called fundamental normal surfaces. 

Algorithms of this type include JSJ decomposition j38] . computing knot genus in arbitrary 
3-manifolds [T| I44j , computing the crosscap number of a knot [12] , and determining which 
Dehn fillings of a knot complement are Haken j33j . 

(iii) For some algorithms, the list .5f is much larger again: typically one must first enumerate all 
fundamental normal surfaces, and then (possibly at great expense) extend this list further 
in some way. Algorithms of this type include finding an essential planar surface in a link 
manifold [32] . or computing the Heegaard genus of a 3-manifold 37 . 

Not all integer points in correspond to normal surfaces: there are additional quadrilateral 
constraints that such points must satisfy. As a result, normal surfaces are only obtained from 
a small subset of the extremal rays or the Hilbert basis of This is a powerful means for 
optimisation that underpins much of this paper. 

Case (i) above is a "best case scenario" for normal surface theory. Vertex normal surfaces 
can be enumerated using highly optimised polytope vertex enumeration techniques, explicitly 
tailored for the setting of normal surface theory [7j [13] . There are practical software implemen- 
tations available [IH 119] , which in turn have enabled new theoretical and experimental advances 

Case (ii) is more problematic. Because of the high dimensionality of the problem and the po- 
tential for super-exponentially many solutions [23] , even state-of-the-art Hilbert basis algorithms 
[H [3] are impractically slow for all but the simplest inputs. Although some special cases have 
been analysed by hand [M] [551 [23 , there are no enumeration algorithms or software implemen- 
tations that exploit the rich structure of normal surface theory, and at present the enumeration 
of fundamental normal surfaces is still considered impractically difhcult. 

It is these issues that we address in this paper. We develop new algorithms to enumerate 
fundamental normal surfaces by merging the structure and constraints of normal surface theory 
with classical Hilbert basis algorithms. We also develop and compare software implementations 
of these algorithms, and show that the best of these is practical for "real" non-trivial inputs by 
using it to compute 164 previously unknown crosscap numbers in the Knotlnfo database [ISjl^l 

^The recent paper [12 also computes new crosscap numbers via normal surfaces, but uses integer programming 
techniques instead. The methods here and in |12| are complementary: each can solve problems that the other 
cannot. 
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Case (iii) above is even more difficult, and many such problems remain severely intractable. 
However, we note that they also benefit from our new algorithms, since the first stage is typically 
to enumerate all fundamental normal surfaces before extending them to a much larger list . 

Hilbert basis enumeration is difficult: it generalises polytope vertex enumeration, which is 
itself NP-hard [23 I3i]- Nevertheless, several families of algorithms exist, as well as efficient 
implementations such as the state-of-the-art software package Normaliz [H [3]- It is widely 
noted that no one algorithm for Hilbert basis enumeration is superior in all settings [H [T71 30] , 
and for this reason we work with three well-known algorithms from the literature: 

• We begin with the primal algorithm described by Bruns, Ichim and Koch [51 [3] and im- 
plemented in Normaliz. In essence, this triangulates the cone into simplicial subcones, 
constructs parallelotopes at the tips of these subcones, locates all integer points within 
these parallelotopes, and then reduces this list to a final basis. 

In Section[3]we describe a primal algorithm for enumerating fundamental normal surfaces. 
The key idea is to enumerate all vertex normal surfaces, and use these to build up the 
small region of ^ that satisfies the quadrilateral constraints. We then run each "maximal 
admissible face" of this region through the primal algorithm of Normaliz, and combine 
the results. 

• We follow with the dual algorithm of Bruns and Ichim [2] , based on work by Pottier [41] 
and also implemented in Normaliz. If we express 'rf in the form {x e M^" | x > 0, Ax = 0}, 
this algorithm inductively builds Hilbert bases for cones defined by the first i rows of the 
matrix A. This generalises the double description method for enumerating extremal rays 
[2T1[39]. 

In Section [4] we describe a dual algorithm for enumerating fundamental normal surfaces, 
by filtering out intermediate points that break the quadrilateral constraints. This in turn 
generalises the filtering technique of Letscher for enumerating vertex normal surfaces [7] . 

• We finish with the Contejean-Devie algorithm 17J, which "grows" potential basis vectors 
by incrementing their coordinates one at a time. Although this is inefficient for high- 
dimensional problems, it has an extremely small memory footprint and so we include it in 
our list. 

In Section[5]we present a variant of the Contejean-Devie algorithm for fundamental normal 
surfaces, using a filtering technique similar to that used in the dual algorithm above. 

Through experimental testing, we find in Section [6] that both the primal and dual algorithms 
for fundamental normal surfaces run significantly faster than generic state-of-the-art Hilbert ba- 
sis algorithms. Of these, the primal algorithm is found to have the most consistent performance. 

Regarding bounds on running times: For the simpler problem of enumerating vertex normal 
surfaces, the best theoretical bounds are still severely exaggerated when compared with real 
behaviour [13] . and we expect the same here. We do not derive explicit bounds in this short 
paper, but we do note that the best known bound on the output size is exp(C'(n^)) [23j, and so 
any bounds on running time should be at least this large. 

All implementations are written using the freely-available software package Regina [U [11] , 
and can be obtained from Regina^s online source repository or the forthcoming version 4.91. 
The primal algorithm also incorporates portions of Normaliz, as described above. 
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2 Preliminaries 



Here we give a very brief overview of triangulations and normal surfaces. In keeping with the 
focus of this paper, we concentrate on the algebraic formulation at the expense of geometric 
insight; for further information the reader is referred to the excellent summary in |23j . 

Throughout this paper, the input for a topological problem is a 3-manifold triangulation 
T, built from n tetrahedra by affinely identifying (or "gluing together") some or all of the 4n 
tetrahedron faces in pairs so that the resulting topological space is a 3-manifold (possibly with 
boundary). 

Such triangulations are more general than simplicial complexes: as a result of the face gluings, 
different edges of the same tetrahedron might be identified together, and likewise with vertices. 
Two faces of the same tetrahedron may even be glued together. This general definition allows us 
to express rich topological structures using few tetrahedra, which is important for computation. 

A normal surface in the 3-manifold triangulation 7~ is a properly embedded surface in 7" 
that meets each tetrahedron in a (possibly empty) disjoint union of normal discs, each of which 
is a curvilinear triangle or quadrilateral. Each triangle separates one vertex of the tetrahedron 
from the others, and each quadrilateral separates two vertices from the others, as illustrated in 



Figure 1(a) Normal surfaces may be disconnected or empty. 





(a) Examples of normal discs 



(b) The seven disc types in a tetrahedron 



Figure 1: Normal triangles and quadrilaterals 



There are seven types of disc in each tetrahedron, defined by which edges of the tetrahedron a 
normal disc meets. These include four triangle types and three quadrilateral types, all illustrated 
in Figure 1(b) The vector representation of a normal surface S' is a vector v(S') G Z^", where 
the 7n elements of v(S') count the number of discs of each type in each tetrahedron. Given 
v(S'), it is simple to reconstruct the surface S (up to an isotopy that preserves the simplices of 
T). 

For each normal surface S, v(S') satisfies the matching equations. These are 3/ linear ho- 
mogeneous equations derived from T, where / is the number of non-boundary faces of T. We 
express these equations as A'v{S) — 0, where ^ is a 3/ x 7n matching matrix. In essence, these 
equations ensure that triangles and quadrilaterals in adjacent tetrahedra can be joined together. 
The matching matrix A is sparse with small integer entries. The normal surface solution cone 
^ is the rational cone = {x e K.^" | Ax = 0, x > 0}, which is a pointed cone with vertex at 
the origin. 

No normal surface S can have two different types of quadrilateral in the same tetrahedron, 
since these would necessarily intersect (contradicting the requirement that S be properly embed- 
ded) . We say that a vector x G M^" satisfies the quadrilateral constraints if, for each tetrahedron 
A of T, at most one of the three coordinates of x that counts quadrilaterals in A is non-zero. 
Running through all n tetrahedra, this gives us n distinct constraints of the form "at most one 



of: 



Xk is non-zero" . The quadrilateral constraints are non- linear, and their solution set is 



non-convex. 
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A point called admissible if it lies in the cone and satisfies the quadrilateral 

constraints. By a theorem of Haken (55], an integer vector x G Z^" represents a normal surface 
if and only if x is admissible. The admissible region of is the set of all admissible points in 

Let ^ C be any pointed rational cone with vertex at the origin, and let ^ n Z'' denote 
aU integer points in ^. Then the Hilbert basis of J^nZ'^, denoted Hilb(^nZ'^) or just Hilb(^) 
for convenience, is a minimal set of integer points that generates all of fl Z'' under addition. 
Equivalently, Hilb(^) is the set of ah x e ^nZ'^ for which, if x = y+z with y, z e ^nZ'', then 
either y = or z = 0. Hilbert bases are finite and unique [25l |48], and have many applications 
[16] ■ They can be defined more generally, but for this paper the setting given here is sufficient. 

A fundamental normal surface is a normal surface S for which, if v{S) = v(T) + v([/) for 
normal surfaces T and U, then either v(T) = or v{U) = 0. The fundamental normal surfaces 
are represented by the admissible points in Hilb('^#'), i.e., those points of Hilb('^) that satisfy 
the quadrilateral constraints. 

A vertex normal surface is a normal surface S for which v(5') lies on an extremal ray of "^lo 
and the elements of v(S') have no common factor^ It is clear that every vertex normal surface 
is also a fundamental normal surface, but in general the converse is not true. 



3 The primal algorithm 



The primal algorithm of Bruns, Ichim and Koch comes in several forms here we describe 

the variant most relevant to us. Let ^ C he a pointed rational cone with vertex at the 
origin. The Bruns-Ichim-Koch algorithm takes as input the extremal rays of J^, and computes 
the Hilbert basis Hilb(,iF). The following is an overview of the procedure; full details can be 
found in |2]. 

1. Use a variant of Fourier-Motzkin elimination to compute the support hyperplanes for 
^ and triangulate ^ into simplicial subcones (cones whose extremal rays are linearly 
independent). 

2. For each simplicial subcone .y, build the semi-open parallelotope spanned by the small- 
est integer vector on each extremal ray of S^, as shown in Figure 2(a) Specifically, if 
the extremal rays are defined by integer vectors Vi , . . . , , we build the parallelotope 
{EA^v. |0< A, < 1}. 





(a) Building a semi-open parallelotope 



(b) All integer points in the parallelotope 



Figure 2: Building parallelotopes in the primal algorithm 



^There are several definitions of vertex normal surface in the literature; others allow common factors [5ll31|. 
or use the double cover if S is one-sided in T 1141 1461 . Our definition here follows that of Jaco and Oertel 1301 . 
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3. Identify all integer points within each semi-open parallelotope, as illustrated in Figure 2(b) 
and combine these into a single list, which generates all of !^ (^ifi. Reduce both the 
intermediate and final lists by removing redundant generators (i.e., vectors that can be 
expressed as non-trivial sums of others in the lists). The final reduced list is the Hilbert 
basis Hilb(^). 

In our case, the pointed rational cone is the normal surface solution cone C R^". Running 
the Bruns-Ichim-Koch algorithm directly over can be slow, since is high-dimensional and 
the triangulation can be extremely large. Moreover, it is wasteful: our aim is to enumerate 
fundamental normal surfaces, and so we do not need all of Hilb('^#'), but just its admissible 
points. 

Our solution, instead of triangulating all of is to just triangulate the admissible region 
of We call a face C ^ admissible if every point x G F is admissible. In particular, we 
note that the vertex normal surfaces correspond precisely to the admissible extremal rays of 
A maximal admissible face of is an admissible face that is not a strict subface of another 
admissible face. 

One can show that the admissible region of is precisely the union of all maximal admissible 
faces [HI [SO]. Note that maximal admissible faces may have different dimensions, and if we 
truncate the origin (which is common to all faces) then the admissible region may even be 
disconnected. 

Our strategy now is to enumerate all maximal admissible faces of '^S' , and use the Bruns- 
Ichim-Koch algorithm to triangulate and compute the Hilbert basis for each. This has the 
advantage that maximal admissible faces are often significantly simpler than the full cone 
they have lower dimension [9 , and are often generated by very few admissible extremal rays [5] . 

The following result underpins our strategy: 

Lemma 1. A vector x G M'^" represents a fundamental normal surface if and only if there is 
some maximal admissible face F C_ '^S' for which x G Hilb(_F). 

Proof. Let x = v(S') where S' is a fundamental normal surface. Since x is an admissible point 
of ^, it must belong to some maximal admissible face F C If x ^ Ililb(F) then there are 
non-zero integer vectors y,z E F for which x = y + z. Since F is an admissible face of 
it follows that y = v(r) and z = v(J7) for some normal surfaces T and U, contradicting the 
fundamentality of S. 

Conversely, suppose that x G Ililb(F) for some maximal admissible face F C Then x 
is an admissible integer vector of and so x = v(S') for some normal surface S. If S is not 
fundamental then v(S') = v(T) -I- v(J7) for normal surfaces T and U with v(T),v([/) 0. 
Since they represent normal surfaces, both v(T),v(?7) G "if, and so any face of the cone 
that contains v(S') must contain both summands v(r) and v(J7) as well. In particular we have 
v(T), v(f7) G F, contradicting the assumption that x = v(T) + v(C/) G Ililb(F). □ 

In our algorithms, we represent faces F C using zero sets [21]. For each face F C'^, the 
zero set of F is a subset of {1, ... , 7n} indicating which coordinate positions are zero throughout 
F. We denote this zero set by Z{F), and define it formally as Z{F) — {k\xk — for all x G F}. 

Because 't^ is of the form {x G K^" | Ax = 0, x > 0}, zero sets effectively encode the support 
hyperplanes for each face: it is then clear that for any two faces F, G C we have Z{F) — Z{G) 
if and only if F = G, and Z{F) C Z{G) if and only if G C F. Zero sets can be stored and 
manipulated efficiently in code using bitmasks of length 7n. 
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The following algorithm makes use of admissible zero sets, which correspond to admissible 
faces of ^€ . We call a set z C {!,..., 7n} admissible if the corresponding zero/non-zero coor- 
dinate pattern satisfies the quadrilateral constraints in K^"; that is, for each tetrahedron A of 
our triangulation, at most one of the three coordinate positions that counts quadrilaterals in A 
does not appear in z. 

Algorithm 2 (Maximal admissible face decomposition). Given the set V of all admissible 
extremal rays of ^ , the following algorithm enumerates all maximal admissible faces of ^ . 

1. Build zero sets for all rays i? G V and store these in the set Si. Initialise an empty output 
set A4, which will eventually contain the zero sets for all maximal admissible faces of . 

2. Inductively build sets S2,S3, . . . as follows. To build the set Sk: 

(a) For each zero set z € Sk-i, find all v ^ Si for which z ^ u and zC\v is admissible. 
For each such v, insert zCiv into Sk- If no such v exists, insert z into the output set 
M. 

(b) Reduce Sk to its maximal elements by set inclusion. That is, remove all z G Sk for 
which there exists some strict .superset z' G Sk. 

(c) If Sk is empty, terminate the algorithm and return the output set A4. 

The key observation is that each set Si, once constructed, contains precisely the zero sets 
for all admissible faces of dimension i. A full proof of correctness and termination will be given 
shortly. 

First, however, we make a further observation regarding zero sets. Recall from standard 
polytope theory that, for any two faces F and G of a polyhedron 3^, the join F V G is the 
unique smallest-dimensional face of ^ that contains both F and G, and that for any face 
H C ^ with F,G CH we have F V G C iJ [49j. 

Lemma 3. For any two faces F, G of the normal surface solution cone , we have Z{F V G) = 
Z{F)nZ{G). 

Proof. Since F y G ^ F,G, any coordinate position that takes a non-zero value in either F or 
G also takes a non-zero value somewhere in F V G. Therefore Z{F V G) C Z{F) D Z{G). 

Let H be the face of obtained by intersecting with the supporting hyperplanes {xi = 0} 
for aU i G Z{F) n Z{G). It is clear that Z{F) n Z(G) C Z{H). Moreover, F,G C H since every 
point in F or G lies on all of the supporting hyperplanes listed above. Therefore F W G C H , 
and so Z{F) D Z{G) C Z{H) C Z{F V G). 

We now have Z{F V G) C Z{F) n Z{G) C Z{F V G), and so Z{F V G) = Z{F) n Z{G). □ 

We can now give a full proof of correctness and termination for Algorithm [2] 

Proof of correctness and termination. We first prove by induction that, once constructed, each 
set Si contains precisely the zero sets for all admissible faces of "^^ of dimension i. 

This is clearly true for i ~ 1, since the 1-dimensional faces of are precisely the extremal 
rays, and we use all admissible extremal rays of to fill Si in step 1 of the algorithm. 

Consider now some i > 1. We assume our inductive hypothesis for both Si-i and iSi, and 
prove it for Si in three stages: 
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(i) Every y ^ Si is the zero set Z{F) for some admissible face _F C of dimension > i. 

Suppose that y is constructed in step 2(a) of the algorithm as y — zDv, where z G iSi-i and 
V G Si. By the inductive hypothesis we have z = Z{G) and v = z{R) for some admissible 
faces G, i? C with dim(G) = i - 1 and dim(i?) = 1. 

By Lemma[3]we have y — Z{G V R). Because step 2(a) of the algorithm only constructs 
admissible zero sets y — z Hv, it follows that G V R must be an admissible face. Because 
step 2(a) only considers pairs for which z ^ v, it follows that R ^ G and therefore G W R 
has strictly higher dimension than G; that is, dim(G V R) > i. 

(ii) For every admissible i- dimensional face F we have Z{F) G Si. 

Let G C F be any [i — l)-dimensional subface of F, and let i? C F be any extremal ray of 
F not contained in G. Since F D G,R we have F D G V R, and since i — 1 = dim(G) < 
dim(G V i?) < dim(F) = i it follows that dini(F) ^ dim(G V R), and hence F = G V R. 
By LemmaHwe then have Z{F) = Z{G) n Z{R). 

Since F is admissible, so are G and R, and it follows from the inductive hypothesis that 
Z{G) G Si-i and Z{R) G Si. Because i? ^ G we have Z{G) ^ Z{R), and we already know 
from F that Z{F) = Z{G) n Z{R) is admissible. Therefore the set Z[F) = Z{G) Ci Z{R) 
is inserted into Si in step 2(a) of the algorithm. 

Furthermore, Z{F) is not removed from Si in step 2(b) of the algorithm. This is because, 
by stage (i) above, any strict superset z' D Z{F) that appears in Si must correspond to a 
strict subface F' <Z F ol dimension > i. This is impossible, since F itself has dimension i. 

(iii) For every admissible j -dimensional face F C with j > i, we have Z{F) ^ Si. 

For any such F, let G C F be any i-dimensional subface of F. Because F is admissible, G 
is also, and by stage (ii) above it follows that Z{G) G Si. Since G is a strict subface of F, 
Z{G) is a strict superset of Z{F), and so even if Z{F) is constructed in step 2(a) of the 
algorithm, it will be subsequently removed in step 2(b). 

Together these three stages establish our inductive claim. 

From here it is simple to see that Algorithm [5] terminates: since the cone has finitely 
many faces, there are only finitely many non-empty sets Sk . Note also that Algorithm [2] does 
not terminate until all non-empty sets Sk have been constructed, since if Sk is empty then there 
are no admissible fc-faces, which means there can be no admissible i-faces for any i > k. 

To finish, we show that the output of Algorithm [2] is correct. Let F C be any maximal 
admissible face of dimension i; by our induction above we have Z{F) C Si. Suppose there is 
some V G Si for which Z{F) ^ v and for which Z[F) H w is admissible. Following the argument 
in stage (i) above, Z{F) n v must be the zero set for some admissible face of dimension >i + l. 
Moreover, this must be a strict superface of F, which contradicts the maximality of F. Therefore 
there can be no such v G Si^ and so we include Z(F) in the output set A4 as required. 

Conversely, consider any output z G M. We must have z G Si for some i, and so 2 = Z{G) for 
some admissible i-dimensional face G G^ . If G is not maximal, there must be some admissible 
(i -|- l)-dimensional superface F Z) G, and some admissible extremal ray i? C F not contained 
in G. Following the argument in stage (ii) above, for the pair Z{G) G Si and Z{R) G Si we 
have Z{G) ^ Z{R) and Z{G) n Z{R) is admissible, which means we would not insert z into the 
output set Ai in step 2(a) of the algorithm. Therefore z = Z{G) for some maximal admissible 
face G C □ 
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Note that each zero set in step 1 can be built by examining any non-zero representative of 
the corresponding ray R. In step 2(a) a set might be constructed from many different {z,v) 
pairs, and so it is important when inserting z Civ into iS^ to avoid duphcates. Each Sk is only 
required for the following inductive step (building Sk+i), and can then be discarded. 

We can now present a full algorithm for enumerating fundamental normal surfaces. 

Algorithm 4 (Primal algorithm for fundamental normal surfaces). To enumerate all funda- 
mental normal surfaces in a given triangulation: 

1. Enumerate all vertex normal surfaces — that is, all admissible extremal rays of 'ia — using 
the algorithms in ^ and which are optimised specifically for normal surface theory. 
Let V be the resulting set of admissible extremal rays. 

2. Using V as input, enumerate all maximal admissible faces of ^S' using Algorithm\^ 

3. For each maximal admissible face F (-^ , identify which rays in V are subfaces of F, and 
pass these rays as input to the Bruns-Ichim-Koch algorithm which will then enumerate 
m\h{F). 

4. To finish, take the union ljHilb(_F) over all maximal admissible faces F. This is pre- 
cisely the set of all vector representations of fundamental normal surfaces in the given 
triangulation. 

In step 3, we can use zero sets to quickly identify which rays are subfaces of F: a ray i? G V 
is a subface of F if and only if Z{R) D Z{F). 

From Lemma [H it is clear that this algorithm is correct. The Bruns-Ichim-Koch algorithm 
does indeed enumerate Ililb(i^) in step 3, because each extremal ray of F is admissible, and so 
the rays in V that are subfaces of F are precisely the extremal rays of F. 

4 The dual algorithm 

We turn now to the dual algorithm of Bruns and Ichim [2], which is based on earlier work of 
Pottier [41] . Again this algorithm comes in several forms, and we describe the variant most 
relevant to us. Consider the cone ^ = {x G M'' I Ax = 0, X > 0}, where A is some rational 
m X d matrix. The dual algorithm takes as input the matrix A, and computes the Hilbert basis 
Hilb(^). 

It constructs a series of cones ^^-^^ . . . , beginning with the non-negative orthant 

= M>Q, and finishing with the target cone ^('") = In general, ^^'^^ is defined as for 
^ above but using only the first k rows of the matrix A. At each stage, the dual algorithm 
uses the previous basis IIilb(^'-'^^^^) to inductively compute the next basis IIilb(,^^'^-'). In this 
sense, the dual algorithm generalises the double description method for enumerating extremal 
rays 

The key inductive step of the dual algorithm is the following [2]. 

Algorithm 5 (Bruns-Ichim-Pottier inductive step). Let 3^ C be a pointed rational cone 
with vertex at the origin, let h G M'' be a rational vector, and define the cones 3^^ = {x G 
^ I h • X > 0}, ^_ = {x G I h • X < 0}, and ^0 = {x G ^ I h • X = 0}. Then, given the input 
basis B — Hilb(,^), the following algorithm computes the Hilbert bases IIilb(^+), Hilb(,^_) 
and Hilb(^o)- 
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1. Initialise candidate bases = {x G _B | h • x > 0} and = {x G i? | h • x < 0}. 

2. Expand these candidate bases: for all pairs x G _B_|_ and y G -B- with h • x > > h • y, 
insert the sum x + y into _B_|_ and/ or _B_ according to whether h • (x + y) > and/ or 
h ■ (x + y) < 0. 

3. Reduce the candidate bases: remove all b G -B+ for which there exists some different 
b' G -B+ with b — b' G ^+ . Reduce i?_ in a similar fashion. 

4- If B^ and i?_ are both unchanged after the most recent round of expansion and reduction 
(steps 2 and 3), terminate with Hilb(^+) = B+, Hilb(^_) = and Hilb(^o) = 

n i?_ . Otherwise return to step 2 for a new round of expansion and reduction. 

In essence, step 2 expands the set of integer points that are generated under addition by 
each individual set _B-|_ and -B„, and step 3 removes redundant vectors from each generating set. 
See [5] for full proofs of correctness and termination, both of which are non-trivial results. 

There are many ways to optimise this algorithm. In practice, one can split -6+ and i?_ into 
three sets according to whether h • x > 0, h • x < or h • x = 0, to avoid duplication along 
the common hyperplanc h • x = 0. The expansion and reduction steps can be partially merged, 
so that we first test each new sum x + y for reduction before inserting it into _B+ and/or B^. 
Several other optimisations are possible: see [21 Remark 16], and the "darwinistic reduction" in 
[1 Section 3]. 

As before, running the dual algorithm of Bruns, Ichim and Pottier directly over the normal 
surface solution cone ^ C K''" can be extremely slow. Like the double description method 
for vertex enumeration, the key problem is combinatorial explosion: even if the final Hilbert 
basis Ililb(^) is small, the intermediate bases Hilb(,^'-'^'') can be extremely large. The high 
dimension of the normal surface solution cone exacerbates this problem. 

Our solution, as with the primal algorithm, is to restrict our attention to just the admissible 
region of This time we do not decompose the admissible region into maximal admissible 
faces; instead we embed the quadrilateral constraints directly into the algorithm itself. In this 
way, a single run through the dual algorithm can compute all admissible Hilbert basis elements 
for 

The idea is simple, and mirrors Letscher's filter for the double description method [7]: when 
expanding candidate bases, ignore sums x + y that do not satisfy the quadrilateral constraints. 

Algorithm 6 (Dual algorithm for fundamental normal surfaces) . To enumerate all fundamental 
normal surfaces in a triangulation with the 3/ x In matching matrix A: 

1. Reorder the rows of A using a good heuristic (as discussed further below). 

2. Initialise the candidate basis B^'^^ with all unit vectors in R^". 

3. Inductively construct candidate bases B^^\ . . . , B^^^"^ as follows. For each i ~ 1, . . . , 3/, 
run a modified Algorithm\^ with input basis B^^~^\ using the vector h G M^" defined by 
the ith row of A. The specific modifications to Algorithm\^ are: 

• In step 2, only consider pairs x, y for which x + y satisfies the quadrilateral con- 
straints. 

• In step 3, instead of testing for b — b' G test for b — b' > and h • (b — b') > 0. 
Instead of testing for b - b' G test /or b - b' > and h • (b - b') < 0. 
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When this modified Algorithm\^ terminates with output and set S*-') = H i3_ . 

4- The final set B^^^^ is then the set of all vector representations of fundamental normal 
surfaces. 

Although we use the Bruns-Ichim-Pottier inductive process in step 3 of our algorithm, we 
do not identify any cone ^ C i?*"" for which the input set i?^*"^) is the Hilbcrt basis. This is 
why we modify step 3 of Algorithm [5] we must remove any reference to the unknown cone 

The order in which we process the rows of A can have a significant effect on performance, by 
affecting the severity of the combinatorial explosion |2ll2l]. In step 1 of Algorithm [6l we reorder 
the rows using position vectors 7^ , which essentially favours rows that begin with long strings of 
zeroes. Position vectors help optimise our quadrilateral constraint filter, and have been found to 
perform well as a sorting heuristic for the related double description method; see [7j for details. 

We can prove correctness and termination for Algorithm |6] as follows. 

Proof of correctness and termination. Let denote the k x 7n submatrix of A obtained by 
taking the first k rows (after the reordering in step 1 of the algorithm), and let ^C^) = {x e 
]^7n = 0, X > 0}. We prove by induction that each candidate basis B^*^^ constructed in 

Algorithm [6] consists of all elements of Hilb(,^'''^'') that satisfy the quadrilateral constraints. 

The initial case is simple: 3^^^^ is just the non-negative orthant K>'^, whose Hilbert basis 
consists of the unit vectors in R*"", all of which satisfy the quadrilateral constraints. This 
matches the initialisation of B'-^'> in step 2 of the algorithm. 

Assume now that i?^*"^) contains all elements of Hilb(,^3^*^'^^'') that satisfy the quadrilateral 
constraints, and consider the construction of i?^*^ in step 3 of Algorithm[6l This involves running 
a modified Algorithm [5] (the Bruns-lchim-Pottier inductive step), with input basis _B(*~^) and 
with the vector h g R^" defined by the ith row of the (sorted) matrix A. 

Let ^+ = {x e ^(^-1) I h • X > 0}, = {x e ^(^-1) I h • X < 0}, and = ^+ n = 
{x g |h • X = 0}. We claim that the modified Algorithm [5] terminates with output 

sets B+ and B^ containing those elements of Hilb(^+) and Hilb(^_) respectively that satisfy 
the quadrilateral constraints. If this is true, then setting i?^*'' B+ D B^ will fill iS^'^ with 
those elements of Hilb(^o) = Hilb(^^*^) that satisfy the quadrilateral constraints, as required. 
This will complete our induction, and the final output B^"^^") will then consist of all admissible 
elements of Hilb(,^^'^-''^) — Hilb('^); that is, all vector representations of fundamental normal 
surfaces. 

It remains to prove our claim above, i.e., that the modified Algorithm [5] terminates with 
output sets B+ and i?_ containing those elements of IIilb(^+) and IIilb(^_) respectively that 
satisfy the quadrilateral constraints. 

Consider running the modified Algorithm [5] and the original Algorithm[5]side-by-side, where 
in the original algorithm we use the full basis Hilb(,f3^(*~^)) as input. Let and B^ denote 
the working sets J5_|_ and B- in our modified algorithm, and let B'^ and J5° denote the working 
sets S+ and in the original algorithm. From Bruns and Ichim P], we know that the original 
algorithm terminates with B° = Hilb(^+) and S° = Hilb(^_). 

We claim that, for as long as both algorithms are running side-by-side, the working sets B™ 
and B^ contain precisely those elements of B"^ and i3° respectively that satisfy the quadrilateral 
constraints. We show this again using induction: 

• It is true when the algorithms begin, since we assume from our "outer induction" on i?^*^ 
that the input basis _B*^*~^) for the modified algorithm contains precisely those elements 
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that satisfy the quadrilateral constraints from the input basis Hilb(^^^* ^-') for the original 
algorithm. 

• We now show that, if at some point and B™ contain the elements of B"^ and B°_ 
that satisfy the quadrilateral constraints, then this remains true after a single round of 
expansion (step 2 of Algorithm [5]) . 

It is clear from our modification in step 3 of Algorithm [5] that we do not insert any new 
sums X + y into B™ or i?™ that break the quadrilateral constraints. Moreover, in the 
original algorithm, any sum x + y that does satisfy the quadrilateral constraints must 
have summands x € B'^ and y G B°_ that satisfy them also (since x, y > 0); therefore 
X e and y G B™, and the sum x + y is also considered in the modified algorithm. 

• Likewise, we can show that if at some point and contain the elements of 5° and 
_B° that satisfy the quadrilateral constraints, then this remains true after a single round 
of reduction (step 3 of Algorithm [5]) . 

First, we note that our modification to the reduction step is sound: for any b,b' G S™, 
we have b — b' G ^+ if and only if b — b' > and h • (b — b') > 0, since we already have 
= ^(*~i)b' = from our outer induction. Next, we observe that a vector b G B^ 
is removed in step 3 of the modified Algorithm [5] if and only if it is removed in the original 
algorithm: this is because b satisfies the quadrilateral constraints, and so any b' for which 
b — b' > must satisfy them also, and therefore b' G -B™ if and only if b' G B'^. Similar 
arguments apply to and B^. 

By induction, and contain precisely those elements of B"^ and 5° that satisfy the 
quadrilateral constraints for as long as both algorithms are running. 

It follows that the modified algorithm terminates no later than the original algorithm, since if 
B"^ and 5° are unchanged under some round of expansion and reduction, then their elements 
and B^ that satisfy the quadrilateral constraints are unchanged also. The modified algorithm 
might terminate earlier, but this will not change the final output: if B^ and i?™ are unchanged 
under some earlier round of expansion and reduction, then running any additional rounds will 
leave them unchanged. 

Either way, we see that the final output sets i3™ and i?™ from the modified Algorithm [5] 
contain precisely those elements of the final output sets B'^ = Hilb(,f3^+) and 5° = Hilb(^_) 
that satisfy the quadrilateral constraints, thus establishing our original claim. □ 

5 The Contejean-Devie algorithm 

Our final algorithm is based on the Contejean-Devie method for enumerating Hilbert bases [17]. 
In essence, we "grow" candidate basis elements out from the origin by incrementally adding 
unit vectors. Like the dual algorithm, the Contejean-Devie method works with a rational cone 
^ = {x G I Ak = 0, X > 0}, taking as input the matrix A and computing the Hilbert basis 
Hilb(^). 

Algorithm 7 (General Contejean-Devie algorithm). To compute the Hilbert basis Hilb(^) as 
described above given the rational matrix A as input: 

1. Initialise an empty output set B, and a "working set" S that contains all unit vectors in 
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2. While S ^ %, update B and S as follows: 

• For all X €z S with Ax. = 0, insert x into B . 

• Create a new working set S' . For all x. Cz S with Ax. =/= 0, if there is no h B 
for which x — b > 0, then insert x + into S' for all unit vectors G M'* with 
Ax. ■ Aei < 0. 

• Replace S with the new working set S' . 

3. The final set B is then the Hilhert basis Hilb(^). 

The algorithm works breadth-first: the fcth iteration of step 2 builds vectors with '^Xi ~ k+1 
by incrementing earlier vectors for which '^Xi = k. The condition Ax ■ Aei < roughly means 
that incrementing the ith coordinate of x will move us towards the null space of the matrix A. 

This breadth-first structure creates some problems. First, there is significant redundancy in 
how vectors are generated: each new vector x' e S" is typically constructed from many different 
source vectors x G S. Moreover, the intermediate working sets S can grow to be extremely 
large. 

Contejean and Devie address these issues by modifying their breadth-first algorithm into a 
depth-first, stack-based algorithm [iTl Section 6]. This stack-based algorithm eliminates redun- 
dancy, and at most d working vectors need to be stored at any one time. 

Although the Contejean-Devie method is no longer state-of-the-art, we include it here be- 
cause it uses different techniques, and because the stack optimisation avoids the severe space 
requirements of the other algorithms: its memory use is essentially dominated by the output 
size only. 

For our setting of normal surfaces, we again embed the quadrilateral constraints directly into 
the algorithm. The modification is simple, and the resulting algorithm, instead of "growing" 
vectors through all of M^" , only grows vectors through the much smaller region of R^" in which 
the quadrilateral constraints are met. 

We present our modified breadth-first algorithm below. The stack-based variant is too 
complex to describe in this short paper (see [17] for details), but the essential modification is 
the same. 

Algorithm 8 (Contejean-Devie algorithm for fundamental normal surfaces). To enumerate all 
fundamental normal surfaces in a triangulation with matching matrix A, run Algorithm as 
above with A as input, but modified so that in step 2 we only consider sums x + that satisfy 
the quadrilateral constraints. The output set B is then the set of all vector representations of 
fundamental normal surfaces. 

Proof of correctness and termination. We use the fact that the general Contejean-Devie algo- 
rithm is known to terminate and correctly compute Hilb(=^3^), as described in Algorithm!?] Recall 
that Algorithm |5] is a simple modification of Algorithm [7] we only consider sums x -I- that 
satisfy the quadrilateral constraints. 

Once again, consider running the modified Algorithm [5] and the original Algorithm [7] side- 
by-side. We denote the working sets B and S in the modified algorithm by B™ and 5'™ , and in 
the original algorithm by B° and 5*°. We claim that, for as long as both algorithms are running 
side-by-side, i?™ and S*™ contain precisely those elements of B° and 5° respectively that satisfy 
the quadrilateral constraints. 

We prove this by a simple induction. The claim is clearly true at the beginning of the 
algorithms, where = B" = % and both 5*™ and S" contain all unit vectors in E^". It also 
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remains true each time we update the sets B and S in step 2 of the algorithm, since we only 
ever consider non- negative vectors: therefore if x + e, satisfies the quadrilateral constraints then 
X must also, and if x — b > for some x that satisfies the quadrilateral constraints then b must 
satisfy them also. 

The modified algorithm terminates no later than the original, since if 5*° = then we must 
have S"^ = 0, but again the modified algorithm might terminate earlier. However, this does not 
affect the output, since if S'™ = at some stage of the modified algorithm then further iterations 
of step 2 will maintain S™" — 0, and will produce no additional output vectors in B™. 

Therefore the final output B"^ of the modified algorithm contains precisely those elements 
of the final output B° — Hilb(,?^) that satisfy the quadrilateral constraints. Since ^ = in 
our case, this gives the vector representations of all fundamental normal surfaces. □ 

6 Performance 

We now compare the performance of the new primal, dual and Contejean-Devie algorithms for 
enumerating fundamental normal surfaces (Algorithms 21 [5] and [5]) for "typical" input triangu- 
lations. 

All algorithms are implemented directly in Regina [H [TT] , except for step 3 of the primal 
algorithm which uses Normaliz [H [3] to run the Bruns-Ichim-Koch algorithm over each maximal 
admissible face. The dual algorithm incorporates the relevant optimisations described by Bruns 
and Ichim ppl , and the Contejean-Devie algorithm uses the optimised stack-based implementa- 
tion |17j . All implementations use arbitrary-precision integer arithmetic. 

As a benchmark, we also compare these new algorithms against a state-of-the-art implemen- 
tation of generic Hilbert basis enumeration. For this we run Normaliz directly over the normal 
surface solution cone to compute the basis Hilb('^), and then identify which of its members are 
admissible. This benchmark comparison shows whether our new algorithms are successful in 
using the context of normal surface theory to their advantage. 



Median running times Worst-case running times 




123456789 10 11 123456789 10 11 

Number of tetrahedra (n) Number of tetrafiedra (n) 



Figure 3: Comparison of running times for different algorithms 
Figure [3] aggregates the running times for all four algorithms over 50 817 distinct input 

^Not all optimisations are found to be effective in the normal surface setting, due to the high dimensionality 
and the quadrilateral constraints. See the full version of this paper for details. 
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triangulations that represent 13400 distinct closed 3-nianifoldsQ The data are grouped by the 
number of tetrahedra n, and within each group we plot the median and worst-case times on a 
log scale. All measurements use single-threaded code running on a 2.93GHz Intel Xeon X5570 
CPU. 

We only plot the Contejean-Devie and benchmark algorithms for n < 4 and n < 5 respec- 
tively, since beyond these limits the algorithms did not finish within a week. Two of the 50 817 
inputs for the dual algorithm also failed to finish in a week (both with n = 11), and these two 
pathological cases are omitted from the dual algorithm plot. 

The results are extremely pleasing: the new primal and dual algorithms run significantly 
faster than the benchmark algorithm, and can process much larger inputs as a result. The median 
times for the primal and dual algorithms are comparable, but the worst-case plot highlights an 
interesting feature: for the dual algorithm, the distribution of running times has a long tail, 
and the worst cases have running times many orders of magnitude slower than the median. In 
contrast, the primal algorithm was consistently fast, with all 50 817 cases finishing in under 
1 second. 

Both the Contejean-Devie and primal algorithms show significantly smaller memory foot- 
prints than the dual and benchmark algorithms. See the full version of this paper for details. 

The success of the primal algorithm may be because the cone typically has only few 
maximal admissible faces and few admissible extremal rays [5], and these are represented by 
integer vectors with small coordinates. As a result, the subproblems passed to the Bruns-Ichim- 
Koch algorithm become relatively simple. Again, we explore these issues further in the full 
version of this paper. 

7 Application to crosscap numbers 

To illustrate the practicality of the new primal algorithm, we use it to compute previously- 
unknown crosscap numbers of knots. The genus of a knot K is an invariant describing the 
smallest-genus orientable surface that the knot bounds; on the other hand, the crosscap number 
C{K) is an invariant describing the smallest- genus non-orientable surface that the knot bounds. 

Crosscap numbers are difficult to compute: of the 2977 non-trivial prime knots with < 12 
crossings, 2661 crosscap numbers are still unknown; in contrast, the genus is known for all of 
these knots [T^ . Although there are specialised methods for certain classes of knots [26l [27j |45] , 
there is no general algorithm yet for computing crosscap numbers. However, we can come close: 
by enumerating fundamental normal surfaces, it is possible to either compute the crosscap 
number of a knot precisely or else reduce it to one of two possible values [Hj . 

Due to space limitations we only give a very brief summary of the computations here. For 
more information on crosscap numbers and how they can be computed, see [12j . The website 
EttpT// www. maths. uq.edu.au/~bab/code/ includes data files giving full details of the com- 
putations below, including all of the relevant triangulations and their vertex and fundamental 
normal surfaces. 

The crosscap number algorithm |12[ Algorithm 4.4] works with triangulated knot comple- 
ments: these are triangulated 3-manifolds with boundary, obtained by removing a small regular 
neighbourhood of a knot from the 3-sphere. The algorithm requires the triangulation T to have 
specific properties (an "efficient suitable triangulation"), which can be tested by enumerating 

*These are the 50 816 minimal orientable triangulations from [H] Table 3], plus the minimal triangulation of 
52 X 51. 
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vertex normal surfaces. It then enumerates all fundamental normal surfaces in T, runs some 
simple tests over each (such as whether the knot bounds the surface, and computing the ori- 
entable or non-orientable genus), and then either determines C{K) precisely or declares that 
C{K) G{t,t+1} for some t. 

We process our knot complements in order by increasing number of tetrahedra, and stop at 
the point where the computations take over two days to run. In total, we run the algorithm 
over 293 knots, whose number of tetrahedra n ranges from 5 to 24 (with mean n ~ 19.7 and 
median n = 21). All but 20 of these computations run in under 12 hours, which is very pleasing 
given that these enumeration problems involve up to 7 x 24 = 168 dimensions. 

The final outcomes are also pleasing: we are able to compute 164 previously-unknown cross- 
cap numbers precisely (for the remaining 129 knots we obtain information that is already known). 
A detailed table of results can be found at http://www. maths. uq.edu.au/~bab/cod.e/ 

It should be noted that the paper [H] also computes new crosscap numbers, but using integer 
programming methods instead. Our techniques here require knots whose complements do not 
have too many tetrahedra, but do not need any prior information about the knots themselves. In 
contrast, the integer programming techniques of |12| can work with very large knot complements, 
but they require knots for which good lower bounds on C{K) are already known. In this sense, 
the two methods complement each other well. 
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